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Abstract 

Prosthetists  today  widely  practice  manual  socket  fitting,  which  produces  subjective, 
inconsistent  results.  To  address  this  problem,  the  Computerized  Anthropometry  Research  and 
Design  (CARD)  Laboratory  is  developing  a  computer-aided  socket  design  system  that  acquires 
ultrasound  datasets  of  an  amputee's  residual  limb,  creates  a  3D  model,  and  helps  identify  load- 
bearing  and  pressure-relief  areas.  This  research  project  focuses  on  providing  3D  visualization  rf 
a  residual  limb  to  support  the  CARD  Laboratory's  efforts.  Creating  the  3D  model  of  the  skin 
and  two  bone  contours  requires  two  major  steps:  segmentation  to  identify  the  objects  of  interest 
and  a  surface  tracking  algorithm  to  generate  the  polygonal  database  of  the  surface  contours. 
Low-level  noise,  incomplete  boundaries,  and  widely  varying  intensities  within  the  images 
present  a  difficult  challenge  to  segmentation  as  well  as  to  the  construction  of  a  3D  model. 
Therefore,  different  techniques  have  been  explored  to  achieve  accurate  segmentation  and 
realistic  3D  rendering  of  unique  ultrasound  images.  Among  various  segmentation  techniques 
tested,  the  enhanced  2D  multiresolution  Bayesian  efficiently  produces  accurate  segmentation  of 
outer  skin  contour  and  bone  locations  of  the  lower  limb.  The  basic  technique  applies  a  filter- 
and-decimate  approach  coupled  with  an  adaptive  clustering  algorithm  and  is  modified  to  use 
anatomical  characteristics  to  detect  and  eliminate  artifacts. 


EVALUATION  OF  SEGMENTATION 


FOR  BONE  STRUCTURES  IN  3D  RENDERING  OF 
ULTRASOUND  RESIDUAL  LIMB  IMAGES 


1  Introduction 


1.1  Background 

For  various  medical  reasons,  doctors  perform  approximately  60,000  lower-extremity 
amputations  every  year  in  the  United  States.  For  greater  mobility,  amputees  are  usually  fitted 
with  a  prosthesis  comprising  an  artificial  foot,  shaft,  and  hand-crafted  socket  [MORI95]. 

The  pressures  of  wearing  a  prosthetic  device  affect  each  area  of  a  below-the-knee  (BK) 
amputee's  residual  limb  differently.  Because  of  this  uniqueness  in  response,  the  prosthetist 
must  decide  which  areas  to  load  with  weight  and  which  areas  to  relieve  from  pressure 
[ENGS92],  [SMrT95].  With  this  concern  for  comfortably  distributing  weight  and  making 
allowances  for  sensitive  areas,  designing  a  prosthetic  socket  depends  largely  on  the  location  d 
the  bone  structures  in  the  residual  limb.  The  areas  around  the  tibia,  the  largest  of  two  bones  in 
the  lower  leg,  tend  to  bear  most  of  the  weight  and  pressure  [MORI95]. 

With  these  factors  in  mind,  a  prosthetist  must  customize  a  socket  so  that  it  comfortably 
fits  the  contours  of  an  amputee's  residual  limb  while  simultaneously  providing  stable  support. 
Prosthetists  perform  custom  fittings  by  hand-crafting  positive  molds  of  residual  limbs  to 
fabricate  sockets  for  a  prosthesis  [SMIT95].  See  Figure  1  for  a  sketch  of  a  positive  mold  with 
built-up  relief  areas. 
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Figure  1  Positive  Mold  of  Residual  Limb  [SAND86] 


Because  conventional  fitting  methods  are  performed  manually,  results  are  subjective, 
inconsistent,  and  greatly  dependent  on  the  prosthetist  [ENGS92],  [SMIT95].  Since  they  rely  on 
the  prosthetist's  knowledge  and  experience,  these  methods  also  usually  require  multiple  fittings 
to  ensure  proper  fit  of  the  prosthetic  socket.  Recurrent  fittings  make  fabrication  of  the  prosthesis 
very  costly.  An  average  prosthesis  costs  $4,000  and  prices  can  range  as  high  as  $15,000 
[MORI95]. 

To  reduce  fittings  and  increase  comfort,  the  Computerized  Anthropometry  Research  and 
Design  (CARD)  Laboratory  at  Armstrong  Laboratories  is  developing  a  computer-aided  socket 
design  (CASD)  software  system  that  performs  socket  fittings  with  higher  accuracy  than 
conventional  methods,  leading  to  greater  cost  savings  and  comfort  for  amputees  [HOUS92].  The 
CARD  Laboratory  is  also  working  with  Sytronics,  Incorporated,  Wright  State  University,  and 
the  New  York  Veterans'  Administration  (VA)  in  this  endeavor. 


1.2  Problem 


The  CASD  software  will  provide  accurate  information  on  the  skin  contours  and  bone 
locations  of  a  BK  amputee's  residual  limb.  Three-dimensional  (3D)  visualization  will  support 
the  CARD  Laboratory's  ongoing  efforts  to  create  the  computer-aided  socket  design  software, 
which  will  resolve  the  problems  associated  with  conventional  fitting  methods. 

My  research  primarily  investigated  bone  segmentation  techniques  for  3D  visualization 
of  low-noise  ultrasound  image  slices.  To  obtain  lower  limb  information,  the  CARD  Laboratory 
chose  a  specialized  ultrasound  imaging  technique  because  of  its  numerous  advantages  over 
other  types  of  image  acquisition.  The  technique  and  its  advantages  are  discussed  in  Chapter  2. 

Because  of  noise  and  other  problems  associated  with  ultrasound  images,  my  research 
focused  heavily  on  exploring  segmentation  techniques  for  the  unique  ultrasound  images. 
Segmentation  techniques  that  extract  the  locations  of  the  bone  structures  and  the  outer  limb 
contour  are  vital  for  generating  accurate  lower  limb  models  for  the  CASD  system.  My  research 
also  provided  a  means  for  producing  polygonal  surface  databases  of  the  3D  surface  model  of  the 
skin  and  bone  contours  for  use  by  the  CASD  system. 

1.3  Methodology 

Three-dimensional  visualization  of  a  residual  limb  is  composed  of  two  major  steps: 
segmentation  and  3D  rendering.  Segmentation  separates  the  objects  in  a  volume  dataset  that 
we  wish  to  view  [BARI93],  [CLAR95],  [PAVL82],  [SAKA95a],  [SAKA95b].  This  process  is  the 
first  operation  that  needs  to  be  applied  to  ultrasound  image  slices.  Specifically  for  this  research, 
the  concern  is  to  separate  and  identify  the  bone  locations  from  the  rest  of  the  limb.  Then,  the 
CASD  software  can  display  the  bones  in  relation  to  the  limb's  surface  contours.  After 
completion  of  segmentation,  surface  rendering  defines  surfaces  of  ol^ects  of  interest  and 
displays  these  surfaces  using  computer  graphics  techniques  [HERM90].  Surface  rendering 
creates  a  3D  model  of  the  residual  limb  from  the  ultrasound  dataset. 
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Many  segmentation  algorithms  already  exist  for  MRI  and  CT  modalities  [BEZD93], 
[CLAR95].  In  comparison,  only  a  few  exist  for  ultrasound  images  due  to  inherent  problems 
associated  with  them,  such  as  high  noise  and  low  resolution.  Fortunately,  the  ultrasound  data 
used  in  this  research  is  unique  because  they  exhibit  lower  amounts  of  these  characteristics. 
Therefore,  segmentation  algorithms  devised  for  MRI  and  CT  images  are  more  readily 
applicable  for  evaluation  on  ultrasound  images.  Still,  the  presence  of  low  noise  levels, 
discontinuous  interfaces  between  different  tissues,  and  wide  intensity  variations  within  tissue 
types  poses  a  tremendous  challenge  [SAKA95a]. 

Because  of  its  lower  computational  and  memory  requirements,  I  chose  surface  rendering 
over  volume  rendering,  a  method  that  uses  each  pixel's  color  and  opacity  information  to  create 
3D  images  [HERM90].  Using  segmented  data,  a  surface-tracking  algorithm  was  executed  to 
generate  a  surface  database  of  polygons  for  a  lower  limb  model.  To  design  an  accurately  fitting 
prosthetic  socket,  the  CASD  system  uses  this  surface  database  to  recreate  a  3D  model  of  a  BK 
amputee's  residual. 

1.4  Overview 

The  remainder  of  this  thesis  is  composed  of  the  following  chapters.  Chapter  2  provides 
details  on  the  acquisition  technique  for  ultrasound  images  identifying  unique  characteristics 
and  potential  difficulties.  Chapter  3  presents  an  overview  of  segmentation  and  3D  visualization 
techniques  with  an  emphasis  on  surface  rendering  methods.  It  also  discusses  the  role  of  3D 
rendering  in  the  CASD  system.  Chapter  4  provides  the  implementation  details  and  results  of 
various  segmentation  techniques.  Chapter  5  discusses  3D  imaging  implementation  and  results. 
Chapter  6  recaps  with  a  conclusion  and  recommendations  for  future  research. 
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2  The  Ultrasound  Image  Dataset 


2.1  Introduction 

Ultrasound  scanning  technology  has  become  very  popular  in  medicine.  The  modality's 
increased  usage  can  be  attributed  to  several  factors.  It  is  non-invasive,  portable,  efficient,  and 
relatively  inexpensive  compared  to  CT  and  MRI  systems  [NELS93],  [SAKA95a].  Because  of 
these  attributes,  ultrasound  imaging  Wcis  chosen  by  the  CARD  Laboratory  to  gather  information 
about  the  surface  contours  and  internal  anatomical  structures  of  an  amputee's  residual  limb. 

2.2  Ultrasound  Image  Acquisition  Technique 

Dr.  Ping  He  and  his  research  group  from  the  Biomedical  and  Human  Factors 
Engineering  Department  at  Wright  State  University  have  developed  a  computer-based 
ultrasound  scanning  system  specifically  implemented  for  acquiring  the  images  of  a  BK 
amputee's  residual  limb.  The  system  used  for  data  acquisition  consists  of  four  parts:  a  water 
tank,  a  patient  seating  device,  a  portable  ultrasound  B-scanner,  and  a  personal  computer.  To 
prevent  movement  during  scanning,  the  seat  was  designed  to  hold  the  thigh  firmly  and  the 
tank  contains  a  limb  rest  [HE94],  [HE96],  [XUE91]. 

The  ultrasound  images  gathered  by  Dr.  He's  system  are  quite  unique.  They  differ  in 
quality  and  appearance  from  typical  images  gathered  from  other  ultrasound  techniques  such  as 
the  example  in  Figure  2.  Dr.  He's  compound  image  technique  reduces  the  level  of  noise, 
thereby  enhancing  details  within  the  images  [HE94],  [XUE91]. 
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Figure  2  Typical  Ultrasound  Image  of  Kidney 


Although  the  method's  resolution  depends  on  a  number  of  factors.  Dr.  He  notes  that  the 
compounding  method  improves  the  lateral  resolution,  which  in  a  typical  ultrasound  varies  with 
the  scanning  depth.  Therefore,  using  multiple  scans  to  reconstruct  the  ultrasound  image  of  the 
lower  limb  tends  to  improve  resolution,  reduce  noise,  and  emphasize  the  details  of  the  outer 
skin  and  bone  structures.  See  Figure  3  for  a  diagram  of  Dr.  He's  apparatus. 
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Water  Tank 


Figure  3  Diagram  of  Compound  B-scan  Method  (Transducer  Rotates  Around  Tank) 


The  compound  imaging  technique  is  ideally  suited  for  this  particular  problem  domain. 
The  lower  limb's  smaller  circumference  allows  the  ultrasound  transducer  to  easily  move  360 
degrees  to  obtain  the  multiple  scans  need  to  reconstruct  the  cross-sectional  image  of  the  residual 
limb.  Also,  the  limb's  simple  internal  structure  prevents  blockage  of  important  objects  of 
interest  especially  the  tibia,  which  is  the  weight-bearing  bone  and  the  fibula,  where  the  muscles 
attach  [HE94]. 

Figure  4  and  Figure  5  show  the  similarities  between  an  accurate  sketching  of  the 
internal  anatomical  structures  and  an  actual  cross-sectional  ultrasound  image  obtained  using  Dr. 
He  s  method.  The  locations  of  the  tibia  and  the  fibula  given  in  Figure  4  correlate  very  closely  to 
the  black  regions  bordered  by  bright  streaks  in  the  ultrasound  image  in  Figure  5.  The  partial 
bright  streaks  surrounding  the  black  regions  show  the  interface  boundaries  between  the  soft 
tissues  and  bone  structures. 
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Figure  4  Anatomical  Cross-section  of  Lower  Limb  [SAND86] 
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Figure  5  Ultrasound  linage  of  Lower  Limb 


Each  ultrasound  image  took  approximately  12  seconds  to  scan  and  an  entire  dataset 
consisting  of  between  40  and  50  images  of  the  residual  limb  takes  between  13  and  15  minutes  to 
produce  [HE96].  The  appendix  at  the  end  shows  a  representative  set  of  the  20  ultrasound 
images  used  for  testing  the  algorithms  implemented  in  this  research.  Each  ultrasound  image 
measures  180  millimeters  (mm)  x  180  mm  and  has  an  interslice  gap  of  6  mm.  These  dimensions 
give  an  aspect  ratio  of  0.7  x  0.7  x  6.  Each  image  contains  256  x  256  pixels  with  grayscale 
intensity  range  of  0  to  255.  This  dataset  is  a  representative  subset  of  the  original  dataset  of  40 
slices  which  consisted  of  180  mm  x  180  mm  slices  with  3  mm  interslice  gaps. 


2.3  Summary 

Compared  to  CT  and  MRI  modalities,  ultrasound  imaging  is  inexpensive,  efficient, 
portable  and  noninvasive.  Because  of  these  advantages,  the  CARD  laboratory  chose  the 
ultrasound  imaging  technique  to  acquire  cross-sectional  image  slices  of  a  BK  amputee's  residual 
limb.  Dr.  Ping  He's  compound  imaging  technique  provided  the  ultrasound  image  slices.  He 
specially  designed  and  developed  this  technique  for  gathering  image  data  on  BK  amputee's 
residual  limbs. 

Although  Dr.  He's  technique  produces  images  with  reduced  noise  and  higher  details, 
typical  ultrasound  problems  due  to  noise  and  incomplete  boundaries  still  exist.  These  problems 
present  a  challenge  in  processing  and  recreating  an  accurate  model  of  the  amputee's  limb.  To 
generate  a  lower  limb  model,  each  slice  in  the  ultrasound  image  dataset  needs  segmentation  to 
extract  the  locations  of  the  bone  structures. 
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3  Image  Processing,  Segmentation,  and  3D  Rendering  Techniques 


3.1  Introduction 

Three-dimensional  (3D)  rendering  algorithms  are  used  to  create  computer-generated 
graphical  models  of  human  anatomical  structures.  These  models  provide  precise  recreations  of 
structures  from  two-dimensional  medical  images.  Visualization  of  3D  models  help  medical 
specialists  to  diagnose  and  analyze  patients'  medical  problems. 

The  medical  images  are  acquired  through  computerized  tomography  (CT),  magnetic 
resonance  imaging  (MRI),  ultrasound,  and  other  medical  modalities.  Of  these  various 
modalities,  ultrasound  is  one  of  three  modalities  capable  of  providing  the  ability  to  visualize 
both  soft  and  bony  tissue  [STYT90].  This  fact,  along  with  its  low  cost,  portability  and  non- 
invasive  nature,  makes  ultrasound  the  best  candidate  for  obtaining  data  about  an  amputee's 
residual  limb  [NELS93].  However,  the  special  characteristics  inherent  in  ultrasound  images, 
such  as  noise  and  lack  of  complete  boundaries,  must  also  be  dealt  with  to  create  accurate 
representations.  To  address  these  problems,  ultrasound  images  are  frequently  preprocessed 
using  image  processing  techniques  to  enhance  image  quality  leading  to  more  accurate 
segmentation. 

Since  my  research  concentrates  on  accurately  creating  a  3D  lower  limb  model,  the  scope 
of  this  review  will  first  briefly  focus  on  image  processing.  Then,  this  review  will  explore 
segmentation,  a  preprocessing  step  for  identifying  objects  of  interest  within  images,  and  the  two 
categories  of  3D  rendering  —  volume  rendering  and  surface  rendering.  Discussion  of 
algorithms  will  be  brief  and  will  provide  high-level  overviews  rather  than  including  the 
specific  steps  in  each  algorithm. 
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3.2  Image  Processing 

Ultrasound  images  have  varying  amounts  of  noise,  or  speckles,  that  hide  important 
details  desired  by  physicians  for  3D  visualization.  Image  processing  involves  the  modification 
of  an  image  to  produce  a  new  image.  Generally,  image  processing  techniques  are  used  to 
manipulate  and  enhance  images  so  that  they  provide  important  information  used  in 
applications  such  as  object  recognition,  edge  detection,  and  image  segmentation  [CAST96]. 

Filtering,  an  image  processing  techniques,  is  frequently  incorporated  into  algorithms  to 
reduce  noise  levels  and  enhance  image  clarity  [SAKA95a],  [SAKA95b].  Filtering  can  be 
performed  in  one  of  two  domains:  spatial  or  Fourier.  Spatially  filtering  an  image  by  passing 
the  filter  over  each  image  pixel  is  equivalent  to  the  convolution  of  images  in  the  Fourier  domain 
[CAST96]. 

In  their  3D  rendering  work,  Sakas  and  Grimm  used  three  different  types  of  filters  to 
achieve  different  goals.  They  used  filters  to  reduce  noise,  smooth  contours,  and  close  gaps  on 
fetal  ultrasound  images  that  are  characterized  by  high  levels  of  noise  [SAKA95b].  After  image 
enhancement,  segmentation  is  necessary  to  separate  different  objects  into  appropriate  tissue 
classes  [SAKA95a],  [SAKA95b],  [SMIT95]. 

3.3  Segmentation 

Due  to  the  presence  of  noise  and  clutter  in  ultrasound  images,  automatic  segmentation 
for  identification  of  objects,  such  as  the  outer  skin,  tibia  and  fibula,  is  difficult.  Gaps  in  the 
boundaries  and  wide  variations  of  intensities  within  tissue  classes  further  complicate  the  ability 
for  unsupervised  segmentation  [HECK96].  For  instance,  without  compensating  for  noise, 
conventional  techniques  could  result  in  misinterpreting  boundaries  of  random  speckle  spots  as 
surfaces,  and  missing  boundaries  at  contour  gaps  [LIN91].  At  the  same  time,  segmentation 
must  produce  accurate  results  in  order  to  correctly  identify  contours  of  the  bone  structures  within 
the  lower  limb. 
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In  his  tutorial  on  3D  imaging,  Udupa  identifies  two  widely  used  classes  of  segmentation 
algorithms:  boundary-based  and  region-based.  Boundary-based  techniques  provide  contours 
of  the  objects  of  interest  while  region-based  ones  produce  binary  images  with  locations  of 
structures  given  by  voxels  marked  by  ones  and  all  other  being  zeros  [UDUP91]. 

Voxels  are  the  three-dimensional  equivalent  of  pixels  [STYT91].  They  are  small  cubes 
identified  by  the  x,  y,  and  z  coordinates  of  their  centers.  Figure  6  shows  a  diagram  of  a  3D 
dataset  of  ultrasound  images  of  a  limb  and  a  voxel  representation  pulled  out  from  a  image  slice. 
Each  voxel's  height  is  equal  to  the  thickness  of  an  image  slice  and  each  voxel  has  a  value  where 
meaning  is  based  upon  useful  data  such  as  intensity  or  density. 


Figure  6  Voxel  Representation  for  Dataset  of  Image  Slices 


Boundary-based  segmentation,  also  known  as  edge  detection,  relies  on  the  detection  of 
rapid  intensity  changes  within  the  image  to  find  edges  separating  different  tissues.  Most 
boundary-based  techniques  rely  on  computing  the  gradient,  a  mathematical  derivative 
function,  that  detects  rapid  change  since  edges  are  normally  characterized  by  high  frequencies 
which  coincide  with  bright  intensity  streaks  within  ultrasound  images  [CLAR95],  [UDUP91]. 
Laplacian  edge  detection  is  another  technique  used  for  boundary-based  segmentation.  The 
Laplacian,  a  scalar  second-derivative  operator,  produces  zero-crossings  at  edges.  Castleman 
gives  other  edge  operators  in  his  book  on  digital  image  processing  [CAST96]. 

In  region-based  approaches,  segmentation  is  achieved  using  a  set  of  properties  that  characterize 
voxels  within  the  same  tissue  type,  and  distinguish  them  apart  from  others.  Udupa  regards 
thresholding  to  be  the  simplest  method.  It  distinguishes  voxels  based  on  a  specified  value. 

Another  common  region-based  approach  involves  clustering  voxels  by  a  set  of 
properties  associated  with  different  regions  [IJDUP91].  Pappas's  adaptive  clustering  algorithm 
has  been  the  basis  for  several  research  efforts  in  the  segmentation  of  MRI  data.  His  algorithm 
includes  spatial  constraints  and  local  intensity  variations.  The  adaptive  clustering  algorithm 
calls  for  moving  a  sliding  window  over  an  image  to  calculate  the  intensity  means  for  that  given 
area.  Using  the  iterated  conditional  mode  (ICM)  approach  presented  by  Besag  in  [BESA86], 
each  pixel  is  inspected  to  maximize  the  conditional  probability  equation  given  as  the  following: 

PUly)  oc  Qxp-lZ[—^iy^-n^  )^]  +  IV^U)]  (1) 

s  s 

where  fl  is  the  intensity  mean  for  class  a:  at  pixel  s,  is  the  local  noise  variance,  and  y„  is 
s  ^ 

the  observed  intensity  mean. 

Equation  1  is  maximized  based  on  the  prior  classification  of  the  pixel  and  its  local  noise 
variance.  Relying  on  a  "pairwise  interaction  model"  with  a  pixel's  four  neighbors,  equation  1 
obtains  the  Gibbs  potential  by  the  summation  over  all  cliques 
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(2) 


c 

where  the  clique  potential  is 


Vc(^i,  ^^0  =  { 


-p  if  x,=Xj 
+P  Otherwise 


(3) 


and  P  is  the  Gibbs  random  field  (GRF)  weight  initially  set  at  1.0  with  the  Ap  =  0.55.  The 
Pappas  algorithm  uses  the  Gibbs  random  field  to  model  spatial  constraint  during  segmentation. 
The  Gibbs  potential  is  based  on  the  GRF  model  and  is  the  a  priori  probability  of  the 
segmentation.  After  a  specified  threshold  of  pixels  converged,  the  cycle  is  repeated  by 
decreasing  the  window  size  again  by  half  until  it  reaches  a  specified  minimum  window  size 
[PAPP92].  Other  researchers  such  as  Chang,  Ashton,  and  Yan  have  used  the  Pappas  clustering 
algorithm  to  develop  segmentation  techniques  for  3D  MRI  scans  of  the  brain  [ASHT95b], 
[CHAN93],  [YAN94]. 

Mathematical  morphology,  a  group  of  image  processing  techniques,  is  also  able  to 
segment  images.  This  approach  is  composed  of  a  set  of  binary-operators  which  can  be  used  in 
different  combinations.  Morphological  operators  use  a  structure  element  (SE),  which  acts 
similar  to  a  convolution  kernel,  over  the  original  image  [CAST96].  The  SE  can  be  of  any  size 
and  shape  containing  O's  and  I's  and  is  usually  designed  based  on  knowledge  of  the  object  of 
interest.  Thomas  and  Peters  used  several  morphological  operators  in  their  research  to  extract 
femur  length  from  fetal  ultrasound  images.  Their  method  used  morphological  operators  to 
reduce  noise,  fill  in  small  holes  in  the  femur,  and  remove  unnecessary  extensions  from  the 
femur  [THOM91].  Other  research  by  Sakas  and  Walters  also  tried  morphological  operators  to 
reduce  noise  and  fill  in  gaps  on  the  surfaces  of  ultrasound  images.  They  were  investigating 
various  types  of  filtering  to  extract  surfaces  from  3D  ultrasound  data  [SAKA95a]. 
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3.4  Methods  for  3D  Imaging 


As  mentioned  previously,  three-dimensional  medical  models  created  by  computers 
result  from  two  primary  categories  of  3D  imaging  techniques:  volume  rendering  and  surface 
rendering  [UDUP90].  Volume  rendering  makes  no  assumptions  concerning  underlying 
structures  or  the  surface,  but  instead,  estimates  structures  using  a  voxel's  light-reflection 
characteristics,  and  color  information  [IJDUP90].  Surface  rendering  displays  images  with 
surface  estimates  based  on  image  data  containing  information  on  the  structures  of  the  human 
body. 

In  general,  although  volume  rendering  tends  to  produce  higher  quality,  photorealistic 
images,  it  is  computationally  expensive  and  requires  more  memory  than  surface  rendering 
[HERM90].  Therefore,  the  choice  in  deciding  between  the  two  approaches  for  3D  visualization 
of  the  residual  limb  for  the  CASD  system  depended  largely  on  the  needs  of  prosthetists.  Since 
prosthetists  are  not  worried  about  the  photorealism,  my  research  investigated  the  use  of  a  voxel- 
based  surface  rendering  technique  that  not  only  generates  images  quickly  but  accurately.  Each 
rendering  technique  is  described  below,  but  volume  rendering  will  only  be  discussed  briefly 
since  this  research  focused  on  surface  rendering  techniques. 

3.4.1  Volume  Rendering  Techniques 

Volume  rendering  methods  have  been  developed  to  display  data  directly  from  gray¬ 
scale  volume  data.  Two  primary  techniques  for  volume  rendering  are  octree  encoding  and  ray¬ 
tracing.  In  octree  encoding,  a  hierarchical  data  structure  must  be  created  to  represent  the 
volume  data.  This  structure  allows  for  direct  manipulation  of  the  3D  model,  for  example,  slicing 
through  images  to  reveal  inner  tissue  [BARI93]. 

Greater  efficiency,  compared  to  the  octree  method,  makes  the  ray-tracing  technique 
more  widely  used  for  volume  rendering  since  no  preprocessing  is  needed  to  set  up  a 
hierarchical  structure.  Ray-tracing  consists  of  three  steps:  the  geometrical  trcinsformation  to 
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align  the  data  with  the  viewpoint,  spatial  segmentation  to  select  only  the  voxels  in  the  volume 
of  interest,  and  finally  gray-level  segmentation  to  display  objects  of  interest  from  the  volume 
dataset  [BARI93]. 

3.4.2  Surface  Rendering  Techniques 

Newer  surface  rendering  techniques  usually  rely  on  voxel-based  methods  rather  than 
triangulation  algorithms  using  planar  contours  to  produce  3D  representations  and  result  in 
higher  accuracy  and  reliability.  Generally,  these  surface  rendering  techniques  are  composed  of 
three  steps:  segmentation,  extraction  and  tiling  of  the  surface  contour,  and  application  of 
perspective  with  shading  [BARI93]. 

Since  segmentation  has  already  been  discussed,  the  discussion  moves  to  the  next  step  of 
surface  tracking.  A  surface-tracking  algorithm  is  executed  to  extract  and  tile  the  surface  contours 
with  polygons  (usually  triangles)  to  create  a  surface  database.  This  approach  starts  by 
identifying  all  voxels  within  a  dataset  that  have  surfaces  cutting  through  them  to  create  a 
surface  database  [BARI93].  In  the  late  1980's,  Lorensen  and  Cline  presented  the  "marching 
cubes"  method  to  detect  voxels  which  have  surfaces  intersecting  them  and  to  estimate  the 
intersecting  polygons  based  on  surface  topology  [LORE87]. 

The  surfaces  extracted  can  then  be  displayed  using  colors,  shading,  and  illumination  to 
create  realistic  3D  images.  To  increase  realism  and  recognition,  colors  that  approximate  the 
anatomical  structures  are  used  to  display  the  3D  model.  Shading  and  illumination  models  are 
used  to  give  depth  and  orientation  cues.  This  final  step,  incorporating  color,  shading,  and 
illumination,  is  not  covered  in  detail  here  since  our  research  is  limited  to  achieving  accurate 
segmentation  and  surface  definitions  of  a  residual  limb. 

3.5  CASD  Systems  and  Role  of  3D  Rendering 

Due  to  advances  in  technology,  computer-aided  socket  design  (CASD)  and  computer- 
aided  manufacturing  (CAM)  research  for  design  and  fabrication  of  prosthetics  and  orthotics  has 


17 


been  rapidly  progressing  over  the  past  ten  years.  Currently,  several  commercial  systems  are 
available  that  use  CASD  and/or  CAM  technologies  [TOPP90]. 

One  such  system  is  the  CANFIT  system  developed  by  the  University  of  British 
Columbia  and  Shape  Technologies,  Inc.  A  clinical  study  compares  this  system  with  the  manual 
form  of  socket  fitting.  The  study  indicated  favorable  results  for  the  CANFIT  system  but  also 
identified  areas  needing  further  improvement  such  as  3D  visualization  of  the  limb  and  data 
acquisition  techniques  to  show  internal  structures  as  well  as  surface  topology  [TOPP90]. 

To  acquire  a  model  of  the  amputee's  limb,  the  CANFIT  system  used  measurements 
obtained  by  the  prosthetist  to  match  the  amputee's  limb  to  the  closest  fitting  models  in  its 
database  of  nine  models.  The  system  then  took  these  models  matching  a  amputee's  limb  to 
extrapolate  and  generate  a  3D  representation.  Because  of  subjectivity  during  measurements 
and  a  limited  database,  data  acquisition  produced  inaccurate  representations  of  the  amputee's 
limb.  As  a  result,  this  process  required  additional  fittings  to  achieve  acceptable  results.  In 
addition  to  the  previously  mentioned  shortfalls,  the  system  was  too  expensive  to  be  used  on  a 
broad  scale  [TOPP90]. 

Another  noteworthy  project  was  at  Sandia  National  Laboratories  in  New  Mexico. 
Researchers  there  developed  a  system  composed  of  ultrasound  imaging  equipment  to  acquire 
data  and  software  composed  of  both  commercial-off-the-shelf  and  customized  software  packages 
to  produce  3D  renderings  of  the  limb  [MORI95].  Since  the  Sandia  system  used  interactive 
segmentation  algorithms  to  create  a  3D  rendering,  there  is  still  a  need  to  address  improvements 
for  faster,  more  automated  segmentation. 

The  CASD  system  that  Armstrong  Laboratory  has  been  creating  will  exhibit  similarities 
to  both  of  the  systems  mentioned  above.  However,  a  number  of  improvements  will  be 
incorporated  into  this  system.  Like  the  system  at  Sandia,  Armstrong's  system  will  use 
ultrasound  imaging  equipment  to  acquire  sequences  of  image  slices  of  the  lower  limb.  The 
system's  software  will  provide  automated  segmentation  to  depict  the  internal  bone  structures 
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and  surface  topology  of  a  residual  limb.  Then  the  CASD  system  will  be  able  to  accurately 
render  a  3D  model  of  a  limb  for  viewing  and  manipulation.  The  CARD  Laboratory  is  also 
exploring  automatic  landmarking  of  pressure  and  relief  areas  on  the  3D  models. 

3.6  Summary 

Algorithms  leading  to  3D  rendering  of  ultrasound  datasets  need  to  consider  the  special 
characteristics  associated  with  them.  Accurate  segmentation  of  the  ultrasound  images  of  the 
residual  limb  is  very  important  in  providing  a  3D  model  for  the  CASD  system.  Since 
numerous  segmentation  algorithms  already  exist  to  handle  CT  and  MR  images,  the  use  of  these 
algorithms  on  ultrasound  images  needs  investigation.  By  specifically  tailoring  these  techniques, 
the  challenges  presented  by  ultrasound  images  can  be  met.  Using  segmentation  combined 
with  surface  rendering  to  generate  a  realistic  3D  model  of  a  limb  provides  better  understanding 
of  the  surface  topology  and  it's  relationship  to  internal  bone  structures.  This  increases  a 
prosthetist's  ability  to  fabricate  more  comfortable  prosthetic  sockets  in  a  shorter  amount  of  time 
with  fewer  fittings. 
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4  Evaluation  of  Segmentation  Techniques 


4.1  Introduction 

This  chapter  presents  the  implementation  details  and  results  of  the  four  segmentation 
techniques  evaluated  in  this  research.  Each  technique's  approach  was  reviewed  for 
applicability  to  low  noise  ultrasound  images.  Since  the  ultrasound  images  used  in  this  research 
contained  low  noise,  these  images  may  not  need  the  extensive  preprocessing  of  other  ultrasound 
images.  Normally,  because  of  high  noise  and  poor  image  quality,  typical  ultrasound  scans 
require  the  use  of  various  image  processing  techniques  to  distinguish  desired  detail  from 
artifacts.  The  examination  of  each  technique  also  ensured  that  algorithms  did  not  rely  on 
parameters  dependent  on  a  particular  modality  such  as  CT  or  MRI.  In  addition,  selection  was 
based  on  consideration  of  image  types  used  for  experimental  testing.  Algorithms  that  had  been 
used  on  images  similar  to  the  ultrasound  lower  limb  images  were  more  favored.  These 
selection  criteria  enabled  overall  assessment  of  each  technique's  potential  for  accurately 
segmenting  the  ultrasound  images. 

The  techniques  implemented  in  this  chapter  were  tested  using  Dr.  Ping  He's  ultrasound 
dataset.  Figure  7  shows  an  ultrasound  slice  from  his  dataset.  The  two  white  outlines  in  the 
image  are  bounding  boxes  denoting  locations  of  the  bone  structures  in  the  lower  limb.  These 
bone  regions  were  manually  identified  based  on  anatomical  knowledge  and  known  ultrasound 
characteristics  of  bone  structures.  To  aid  in  comparison,  most  figures  in  the  result  sections 
display  composite  images.  A  composite  image  has  the  original  ultrasound  image  overlaid  with 
the  boundary  outlines  from  the  segmented  image. 
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Figure  7  Locations  of  Bone  Structures  in  Original  Image 


4.2  MultiresoluHon  Bayesian  Segmentation 

Examination  of  previous  research  efforts  in  the  use  of  a  Bayesian  approach  provided 
evidence  for  its  potential  success  in  segmentation  of  noisy  images.  In  an  earlier  research  effort, 
the  Bayesian  technique  accurately  detected  and  analyzed  surface  classifications  and  normals  on 
noisy  ultrasound  images  [LIN91].  In  two  other  efforts,  the  Bayesian  approach  provided  accurate, 
automated  segmentation  of  MRI  brain  scans  [CHAN93],  [YAN94].  According  to  Sakas  and 
Walter,  a  multiresolution  framework  allows  for  elimination  of  noise  and  minor  structures  faster 
than  important  surfaces.  This  happens  because  smaller  artifacts  have  lower  coherency  while 
larger  surfaces  are  still  detectable  at  lower  resolutions  [SAKA95a]. 

Supported  by  these  previous  efforts,  the  multiresolution  Bayesian  approach  presented 
by  Ashton  and  Parker  was  selected  for  its  ability  to  reduce  speckles  and  accurately  segment 
noisy  2D  ultrasound  images.  The  developers  of  this  algorithm  tested  their  approach  on 


transesophageal  echocardiographic  scans  of  the  human  heart.  The  specialized  algorithm  takes 
into  account  the  statistical  nature  of  ultrasound  data  to  process  and  segment  images  [ASHT94a], 
[ASHT95a].  The  following  paragraphs  summarize  the  implementation  of  this  algorithm. 

4.2.1  Approach  and  Implementation  Of  Original  Algorithm 

This  Bayesian  algorithm  initially  applies  a  low-pass  filter  to  remove  speckles  from  the 
original  ultrasound  images.  The  filter,  constructed  to  characterize  and  reduce  the  speckles,  has 
to  optimally  maximize  two  goals:  noise  elimination  and  preservation  of  important  edge  detail. 
The  filter  suggested  by  Ashton  and  Parker  is  a  finite  impulse  response  (FIR)  filter.  The  FIR 
filter  is  implemented  as  a  Gaussian  with  a  response  equal  to  twice  the  diameter  of  the  estimated 
speckle  spot  size  n. 

To  increase  speed,  the  FIR  filter  is  implemented  in  the  Fourier  domain.  First,  an  image 
of  the  filter  is  constructed  using  a  multivariate  Gaussian  distribution.  The  original  image  and 
the  filter  image  have  to  be  the  same  image  size  with  dimensions  based  on  powers  of  two. 
Then,  both  images  are  fast  Fourier  transformed  and  convolved.  The  convolved  image  are  then 
inverse  Fourier  transformed  back  to  the  spatial  domain  to  produce  a  filtered  version  of  the 
original  image. 

The  speckle  spot  size  n  is  an  important  filter  parameter.  Varying  the  size  of  the  speckle 
spot  parameter  affects  the  resulting  image.  The  larger  the  speckle  spot  size  (or  diameter  of  the 
Gaussian),  the  lower  the  frequencies  passed.  In  other  words,  a  smoother  image  results. 
Increasing  the  smoothing  factor  of  the  filtered  image  can  cause  the  loss  of  important  boundary  or 
edge  definition.  Therefore,  the  speckle  spot  size  has  to  be  carefully  chosen. 

After  filtering,  the  image  is  segmented  with  a  K-means  clustering  algorithm  to  calculate 
each  pixel's  initial  intensity-based  classification.  The  K-means  algorithm  finds  the  k  number  of 
intensity  means  that  an  image's  intensity  values  cluster  around.  These  k  intensity  means  are 
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used  to  classify  each  pixel  in  the  image.  This  initial  classification  provides  the  basis  for  the 
adaptive  clustering  algorithm. 

Ashton  and  Parker  created  a  modified  version  of  the  Pappas  adaptive  clustering 
algorithm  summarized  in  Chapter  3  to  use  as  the  backbone  of  their  algorithm.  Their  modified 
adaptive  clustering  algorithm  classifies  each  pixel  based  on  local  averaging  using  the  intensity 
means  of  surrounding  pixels.  The  classification  of  the  pixeTs  four  neighbors,  which  form  four- 
pairwise  cliques  with  the  pixel,  is  also  considered  as  well  as  its  local  noise  variance  when 
maximizing  the  conditional  probability  function  given  in  Equation  4. 

The  adaptive  clustering  algorithm  is  initially  executed  on  the  lowest  resolution  image 
resulting  from  the  decimation  of  the  filtered  image  by  the  original  speckle  spot  size  n.  After 
maximization  of  Equation  4,  this  process  continues  to  the  next  higher  resolution  by  replicating  to 
the  original  image  size  followed  by  decimation  where  n  =  n  -  1.  Upon  reaching  the  highest 
resolution,  equivalent  to  the  original  image  size,  a  final  execution  of  the  adaptive  clustering 
algorithm  produces  the  final  pixel  classifications  for  the  segmented  image. 

Equation  1  given  in  Chapter  3  is  modified  to  maximize  the  conditional  probability 
equation  given  as  the  following: 

Pix\y)  oc  n(— )exp-[X[^2-(>’5--W;c  +  (4) 

s  ^ X  s  2(J  ^  c 

where  is  the  intensity  mean  for  class  x  at  pixel  s,  CT  is  the  local  noise  variance,  and  y  is 
s  ^ 

the  observed  intensity  mean.  This  modification  by  Ashton  and  Parker  accounts  for  the  fact  that 

ultrasound  images  are  not  corrupted  by  globally-uniform  Gaussian  noise.  However,  since 

Gaussian  statistics  are  easier  to  represent  mathematically  than  Rayleigh  statistics,  they  assumed 

a  Gaussian  distribution  for  pixels  within  a  given  region  where  (T  ,  the  local  noise  variance,  is 

estimated  independently  for  each  class  and  is  proportional  to  the  local  class  mean.  After  a 
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specified  number  of  pixels  converged,  the  local  noise  variance  is  recalculated  and  the  cycle  is 
repeated  by  decreasing  the  window  size  again  by  half  until  it  reaches  a  specified  minimum 
window  size.  The  general  steps  of  the  algorithm  are  as  follows: 

1.  Low-pass  filter  the  original  image  with  a  Gaussian  filter 

of  diameter  2n  x  2n.  Decimate  image  by  n  (the  size  of  speckles). 

2.  Obtain  initial  cluster  estimates  by  K-means. 

3.  Set  the  GRF  weight  to  an  initial  value. 

4.  Apply  the  clustering  algorithm. 

5.  Expand  the  cluster  matrix  by  n  and  then  decimate  by  n-1. 

6.  Decrease  n  by  1. 

7.  Lowpass  filter  and  decimate  the  original  image  by  new  n. 

8.  If  n  >  1  then  increase  GRF  weight  and  go  to  step  4. 

9.  Apply  clustering  algorithm  a  final  time  on  original  image. 


4.2.2  Results  of  Original  Algorithm 


The  original  2D  multiresolution  Bayesian  approach  produced  a  segmented  image  that 
contained  artifacts  within  the  limb  and  did  not  correctly  distinguish  bone  structures  from  the 
artifacts.  Because  of  noise  and  wide  intensity  variations  throughout  the  images,  segmentation 
identified  numerous  regions  within  the  limb  as  bone  structures.  Therefore,  this  segmentation 
does  not  accurately  identify  bone  locations  which  are  vital  for  generating  3D  surface  models  of 
outer  skin  and  two  bone  contours  in  a  lower  limb.  Figure  8  shows  a  composite  of  an  original 
image  with  contours  of  a  segmented  image. 
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Figure  8  Segmented  Image  Using  Original  2D  Bayesian 


4.2.3  Approach  and  Implementation  of  Enhanced  Algorithm 

Morphological  operators  based  on  shape  characteristics  are  frequently  used  in  image 
processing.  As  described  in  Chapter  3,  a  structure  element  (SE)  acts  similar  to  a  convolution 
kernel  and  is  passed  over  each  pixel  in  the  original  image.  The  morphological  operation  of 
dilation  fills  in  holes  in  the  image  that  are  equal  to  or  smaller  than  the  size  of  the  structure 
element.  The  morphological  operation  of  erosion  removes  small  clusters  that  are  smaller  than  a 
given  structure  element.  When  erosion  follows  dilation  of  images,  this  is  known  as  a  closing 
operation  [LYBA94].  Combined  in  this  meinner,  these  operators  can  help  remove  noise  or 
artifacts,  and  refine  contours  producing  better  segmentation  than  algorithms  that  lack 
morphological  operators  [SAKA95a],  [THOM91].  Due  to  these  benefits,  this  research  developed 
new  enhancements  for  use  in  conjunction  with  the  2D  Bayesian  technique.  Morphological 
operators  were  combined  with  the  2D  Bayesian  to  extract  of  bone  regions  from  artifacts. 

To  enhance  Bayesian  segmentation  of  the  ultrasound  image,  a  binary  image  of  each 
image's  potential  bone  regions  is  created  by  applying  a  mask  to  the  segmented  image.  Then 


execution  of  the  closing  morphological  operator  removes  small  artifacts  and  refines  contours  of 
candidate  bone  regions.  The  remaining  image  of  potential  bone  regions  requires  further 
processing  with  additional  discriminators  to  distinguish  the  tibia  and  fibula  from  artifacts. 
These  discriminators  remove  regions  not  meeting  certain  size,  shape,  and  location  criteria 
associated  with  the  two  bone  (tibia  and  fibula)  regions. 

To  extract  the  two  bone  structures,  the  enhanced  Bayesian  algorithm  uniquely  labels 
each  region  of  contiguous  pixels  by  looking  at  their  8-connected  neighbors.  Using  the  labels, 
this  algorithm  then  calculates  each  region's  size  and  identifies  the  largest  region  as  the  tibia. 
After  the  selection  of  the  tibia,  the  algorithm  searches  for  potential  fibula  candidates  by 
checking  the  remaining  candidates'  sizes  to  see  if  they  are  at  least  1/6  the  size  of  the  tibia.  This 
ratio  was  chosen  after  evaluating  several  values. 

For  further  evaluation  of  potential  fibula  candidates,  the  algorithm  compares  each 
candidate's  distance  from  the  tibia  for  closeness  and  shape  for  circularity.  It  calculates  the 
distance  from  the  center  points  of  the  bounding  boxes  for  the  tibia  and  fibula  regions  and 
weights  this  by  0.7.  Circularity  is  assigned  a  weight  of  0.3.  These  weights  were  determined  by 
observation  of  the  ultrasound  dataset  and  may  need  adjustments  with  further  testing  of  more 
ultrasound  datasets.  Examining  initially  segmented  images,  the  minimum  distance  of  the 
fibula  region  to  the  tibia  is  more  important  than  the  region's  overall  roundness.  The  final 
candidate  having  the  minimum  weighted  distance  and  shape  is  identified  as  the  fibula. 

Following  identification  of  bone  structures,  the  enhanced  implementation  assigns  pixels 
of  each  tissue  class  a  specific  intensity  of  50  for  the  limb  and  255  for  each  of  the  two  bones. 
These  values  are  used  later  in  chapter  5  as  the  isovalues  for  the  "marching  cubes"  algorithm. 

For  fast  prototyping  and  quick  results,  the  enhanced  multiresolution  Bayesian  approach 
was  implemented  using  MATLAB  version  4.2  on  the  SGI  workstations.  MATLAB  code  was 
written  to  input  segmented  images  and  perform  the  extraction  of  bone  regions.  The  code  used 
built-in  morphological  operators  to  perform  the  closing  operation  with  a  structure  element  10 
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voxels  in  diameter.  All  enhancements  except  for  the  closing  operation  were  incorporated  into 
the  original  Bayesian  segmentation  algorithm  in  C++. 

4.2.4  Results  of  Enhanced  Algorithm 

The  modified,  enhanced  version  evaluated  and  distinguished  the  tibia  quite  accurately  . 
The  results  were  verified  by  Dr.  Ping  He  and  Dr.  Kefu  Xue,  two  biomedical  experts  in 
ultrasound  imaging  from  Wright  State  University.  They  also  observed  that  segmentation  of  the 
tibia  appeared  quite  accurate.  The  tibia's  accuracy  is  important  for  an  amputee's  residual  limb 
since  this  bone  is  load-bearing. 

Without  the  closing  morphological  operator,  the  tibia's  back  shadow  in  the  ultrasound 
image  detracted  from  segmentation  of  the  fibula  and  led  to  elongated  and  misshapen  contours. 
Figure  9  displays  a  segmented  image  without  a  closing  operator.  Inclusion  of  the  closing 
morphological  operator  in  the  enhanced  Bayesian  algorithm  helped  to  refine  the  contours  of  the 
bones.  In  comparison  to  Figure  9,  the  composite  image  in  Figure  10  shows  more  accurately 
defined  bone  contours  resulting  from  the  closing  operator.  Appendix  B  and  Appendix  C  show 
the  entire  segmented  dataset  of  ultrasound  images  using  the  original  and  enhanced  algorithms. 


Figure  9  Enhanced  Segmented  Image  Without  Cbsing 


Figure  10  Enhanced  Segmented  Image  With  Closing 


4.3  Three  Dimensional  Bayesian  Segmentation 

Three  dimensional  (3D)  segmentation  takes  into  account  the  relationship  of  a  pixel  to  its 
adjacent  slice  neighbors  as  well  as  its  intraslice  neighbors.  As  a  result,  3D  segmentation  allows 
connectivity  and  smoothness  constraints  to  be  imposed  in  3D  and  may  increase  accuracy 
[CHAN93].  Therefore,  by  modifying  the  original  algorithm  described  in  section  4.2  for 
multiresolution  Bayesian  to  include  3D  segmentation,  the  ultrasound  dataset  may  be  more 
accurately  segmented. 


4.3.1  Approach  and  Implementation 

The  3D  Bayesian  technique  is  similar  to  the  previously  described  2D  version.  The  3D 
Bayesian  technique  relies  on  the  modified  adaptive  clustering  algorithm  and  Besag's  ICM  to 
maximize  a  pixel's  conditional  probability  based  on  its  relationship  with  all  neighboring  pixels 
between  and  within  slices  [ASHT94b],  [ASHT95b],  [CHAN93],  [YAN94].  Unlike  the  previous 
technique,  the  3D  Bayesian  approach  does  not  apply  filtering  before  segmenting  each  image. 
However,  to  take  into  account  the  Gibbs  potential  between  neighboring  interslice  pixels,  the 
clique  potential  is  extended  to  3D  as  follows: 


I  A 

if  X.=Xj 

otherwise 

(5) 

f-A 

ii 

(6) 

otherwise 

with  the  Gibbs  potential  given  by  the  sum  of  all  clique  potentials  between  the  four  intraslice 
neighbors  and  two  adjacent  slice  neighbors.  The  intraslice  Gibbs  weight  is  set  at  Pi  =2.0  and  the 
adjacent  slice  Gibbs  weight  at  P2  =0.12.  These  weights  are  calculated  based  on  an  intraslice 
resolution  of  0.7  mm  and  an  adjacent  slice  resolution  of  6.0  mm  using  the  ratios  given  in 
[CHAN93].  Therefore,  the  pixels  are  more  likely  to  be  in  the  same  class  as  their  intraslice 
neighbors  than  their  interslice  ones. 
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Like  the  2D  version,  the  3D  approach  initially  classifies  the  pixels  in  each  slice  by 
applying  K-means  clustering.  Then  for  each  slice  of  the  3D  dataset,  the  adaptive  clustering 
algorithm  is  maximized  after  substituting  the  extended  clique  potentials  given  by  Equation  4 
and  Equation  5  into  the  conditional  probability  formula  from  Equation  1  [CHAN93]. 

43.2  Results  of  3D  Bayesian 

Applying  the  3D  Bayesian  approach  did  not  yield  better  results  than  the  2D  method 
because  of  the  low  interslice  resolution  between  slices  and  the  presence  of  noise  in  the 
ultrasound  images.  The  image  in  Figure  11  shows  the  segmented  image  produced  by  the 
algorithm  and  shows  the  increased  presence  of  artifacts  due  to  noise. 

Preprocessing  all  ultrasound  slices  in  the  dataset  with  a  FIR  filter  (a  Gaussian  with 
diameter  of  2n  x  2n  where  the  speckle  spot  size  is  n  =  2)  produced  segmented  images  that 
closely  resembled  images  from  the  original  2D  Bayesian.  The  filter  reduced  the  level  of  noise 
in  the  ultrasound  slices  and  reduced  the  artifacts  present  in  the  segmented  data.  Figure  12 
shows  a  filtered  image  from  the  3D  dataset  that  looks  similar  to  the  image  slice  1  given  in 
Appendix  B. 
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Figure  11  3D  Bayesian  Segmented  Image 


Figure  12  Segmented  Image  Using  3D  Bayesian  with  Filtering 
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4.4  Marr-Hildreth  and  Morphological  Operators 

Edge  detection  algorithms  can  be  applied  during  segmentation  of  MR  images.  Bomans' 
work,  which  combines  edge  detection  and  morphological  operators,  was  chosen  to  explore  the 
effectiveness  of  boundary-based  techniques  on  ultrasound  images.  Bomans  used  a  3D  Marr- 
Hildreth  operator  along  with  the  dilation  and  erosion  morphological  operators  to  segment 
magnetic  resonance  images  of  a  human  head  for  3D  rendering.  As  mentioned  during  the 
discussion  of  the  enhanced  Bayesian  technique,  the  use  of  morphological  operators  tends  to 
refine  contours  while  removing  artifacts.  Bomans'  technique  resulted  in  accurate  segmentation 
of  MRI  brain  [BOMA90]. 

4.4.1  Approach  and  Implementation 

Using  zero  crossings,  the  locations  of  intensity  changes,  the  Marr-Hildreth  edge  operator 
can  detect  anatomical  structures  contained  in  images.  The  Marr-Hildreth  operator  is  defined  by 
the  following: 

E{x,y,z)  =  V^(/(x,y,z)*G(x,y,z,cr))  (7) 

where  V  is  the  Laplacian  operator,  l(x,y,z)  is  the  image,  G(x,y,z,a)  is  the  Gaussian,  and 
E(x,y,z)  is  the  image  of  edges.  The  operator  smoothes  the  data  to  remove  areas  of  high 
frequency  and  uses  a  Laplacian  to  detect  edges  regardless  of  orientation  [CAST96].  Then  3D 
morphological  operators  of  dilation  and  erosion  are  applied  to  enhance  surface  contour 
definition.  With  these  surface  contours,  Bomans  generated  a  realistic  3D  display  of  the  head 
consisting  of  skin,  bone,  brain,  and  the  ventricles  [BOMA90].  This  implementation  does  not  go 
into  the  other  steps  used  by  Bomans,  et  al.  to  obtain  the  3D  rendering  since  the  area  of  focus 
was  the  comparison  of  segmentation  approaches  on  this  type  of  ultrasound  images. 

MATLAB  was  used  to  generate  the  code  and  test  this  approach  on  the  ultrasound 
images  of  the  lower  limb.  Slight  modification  of  MATLAB's  built-in  functions  of  the  2D  Marr- 
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Hildreth  edge  operator  and  the  2D  closing  morphological  allowed  for  the  testing  of  their 
effectiveness  in  segmenting  the  ultrasound  images.  The  Marr-Hildreth  operator  was  modified 
to  approximate  a  Laplacian  of  a  Gaussian  using  the  Difference  of  Gaussians  as  given  below 
[BOMA90].  Also,  the  MATLAB's  closing  operator's  structure  element  was  modified  to  a  disk 
shape  of  diameter  5.  Although  Bomans  extended  these  operators  to  three  dimensions,  the  two 
dimensional  operators  in  MATLAB  were  able  to  provide  a  general  idea  of  how  this  approach 
would  segment  ultrasound  images.  If  the  2D  methods  in  MATLAB  provided  better 
segmentation  of  the  skin  and  bone  regions  than  the  previous  techniques,  then  implementation 
in  3D  may  yield  similar  or  better  results. 

This  implementation  is  similar  to  Bomans  approach  of  using  a  difference  of  Gaussian 
(DOG)  operator  to  estimate  the  Laplacian  of  a  Gaussian  or  V  G  operator  defined  as 

V^G(x,  y,  z,  (t)  =  G{x,  y,  z,  <Tg )  -  G{x,  y,  z,  ct,  )  (8) 

Bomans  used  the  DOG  operator  to  reduce  memory  requirements  and  lower  the  number  of 

convolutions.  If  a  DOG  operator  is  used  to  approximate  the  V  G  operator  then  the  standard 
deviations  and  <7^  are  set  to  make  the  filter  bandwidth  small  and  the  sensitivity  large. 
The  ratio  of  was  1.6  where 


"  V  21nr 

(9) 

(10) 

and  r  -  1.6  and  —  0.805 -CT  and  CTi  —1.288*0'  [BOMA90].  This  implementation 

attempted  standard  deviations  of  (T  =  1.0,  2.0,  and  3.0  to  examine  the  difference  in  the  resulting 
images.  The  best  results  came  from  using  standard  deviation  of  1.0.  For  the  morphological 
operators,  dilation  and  erosion,  a  structuring  element  of  both  5-voxel  diameter  and  10-voxel 


33 


diameter  was  applied  to  close  the  binary  images  from  the  Marr-Hildreth  operator  to  eliminate 
small  holes  and  smooth  out  contours. 

4.4.2  Results  of  Marr-Hildreth  and  Morphological  Operators 

Using  the  algorithm  developed  by  Bomans  with  a  structure  element  (SE)  of  diameter  5 
did  not  produce  good  results.  The  structure  element  was  not  big  enough  to  eliminate  small 
artifacts  after  the  edge  operation.  As  a  result,  the  final  image  had  artifacts  and  the  technique 
did  not  detect  the  bone  structures  as  well  as  the  enhanced  2D  Bayesian  approach.  Figure  13 
shows  two  images:  the  top  shows  all  edges  detected  using  the  Marr-Hildreth  operator  and  the 
bottom  shows  the  new  image  cifter  the  closing  operation.  Since  the  morphological  operators  rely 
on  size  and  shape  characteristics,  the  SE  diameter  was  increased  to  10  to  remove  the  small 
circular  artifacts  shown  by  the  lower  image  in  Figure  13.  This  size  eliminated  most  of  the 
artifacts  and  left  two  regions  remaining  that  correspond  closer  to  the  locations  of  the  bone 
structures.  The  final  segmented  image  was  similar  to  the  one  obtained  from  the  enhanced  2D 
Bayesian  technique.  However,  this  approach  introduced  two  small  protruding  artifacts  along 
the  upper  skin  contour  of  the  limb  due  to  the  presence  of  noise  clusters.  Figure  14  shows  the 
result  of  reapplying  the  edge  operator  and  the  closing  operator  with  SE  of  10.  The  enhanced  2D 
Bayesian  technique  provides  better  segmentation  without  introducing  additional  artifacts  in  the 
final  segmented  image. 
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Figure  13  Marr-Hildreth  and  Closing  with  SE  5  Figure  14  Marr-Hildreth  and  Closing  with  SE  10 


35 


4.5  Contour  Searching  with  Genetic  Algorithm  (G A) 

Heckman  and  Krumm  of  Sandia  National  Laboratories  developed  a  boundary-based 
segmentation  algorithm  designed  for  ultrasound  images  of  the  lower  limb  (similar  to  Dr.  He's 
images).  According  to  their  research,  automatic  segmentation  of  these  images  is  difficult 
because  of  various  amounts  of  noise  and  "complicated  interactions  between  ultrasonic  pulses 
and  human  tissue"  [HECK96]. 

Edge  finding  and  region  growing  techniques  do  not  work  well  and  without 
postprocessing  tend  to  produce  artifacts  caused  by  noise  and  incomplete  boundaries.  Heckman 
and  Krumm  also  state  that  snakes,  a  contouring  method  based  on  spatial  derivatives,  did  not 
work  well  with  ultrasound  images  and  required  an  accurate  starting  position.  The  presence  of 
noise  diverted  the  snakes'  contours  from  the  actual  boundaries  in  the  images  [HECK96]. 

To  overcome  these  difficulties,  the  authors  used  a  genetic  algorithm  (GA).  The  GA 
searches  for  the  best-fitting,  closed  contours  given  by  a  closed  cubic  spline  to  the  original  image 
of  the  outer  skin  and  bone  structures.  Their  approach  is  similar  to  the  Hough  transform,  a 
powerful  image  processing  technique  used  to  detect  edges,  because  it  globally  searches  for  the 
correct  contours.  However,  their  technique  is  also  able  to  conform  to  the  complex  shapes  present 
in  the  images  [HECK96]. 

4.5.1  Approach  and  Implementation 

The  GA  tries  to  optimize  a  fitness  function  in  its  search  for  the  best  fitting  contour  in  a 
very  large  search  space.  The  algorithm  represents  the  fitness  variables  with  a  string  of 
symbols.  Then,  the  GA  randomly  creates  a  population,  a  set  of  the  strings,  to  begin  the  search 
and  the  fitness  function  is  calculated  for  each  member  of  the  population.  Mutation,  crossover, 
and  selection  techniques  are  also  applied  to  each  member  until  one  member  obtains  the 
appropriate  fitness  level.  Specifically,  mutation  produces  the  global  optimum;  crossover  creates 
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high  fitness  members  using  parts  of  other  high  fitness  strings;  and  selection  tries  to  choose  high 
fitness  strings  by  using  a  biased  probability  distribution  to  pick  them  randomly  [HECK96]. 

The  closed  cubic  splines  used  to  model  contours  have  four  control  points  which  allow  for 
fairly  smooth  contours  while  reducing  the  search  space.  Since  the  search  space  is  still  quite 
large,  Heck  and  Krumm's  algorithm  aligns  initial  control  points  to  likely  boundaries  or  high 
intensity  regions  in  the  image  to  speed  up  the  search  process.  Before  correlating  the  contours  to 
the  original  image,  each  contour  is  refined  to  fall  on  "peaks  of  nearby  intensity  ridges"  of  the 
original  data  to  compensate  for  discontinuities.  This  refinement  process  places  perpendicular 
search  lines  along  the  contours  to  locate  maximum  intensity  values  that  are  used  to  align 
contours.  After  refinement,  contours  are  modified  with  a  point  spread  function  to  simulate  how 
they  would  look  in  a  real  ultrasound  image  [HECK96].  To  decrease  the  search  space  and  make 
the  computing  less  intensive,  the  authors  designed  the  GA  to  search  each  contour  separately.  In 
addition,  the  GA  requires  a  known  bounding  box  for  each  contour's  location  so  that  the 
algorithm  does  not  incorrectly  chose  large  areas  of  low  intensity  with  high  fitness  as  the  best 
contours  [HECK96]. 

The  C  source  code  for  this  algorithm  was  obtained  from  the  authors  and  compiled  on  a 
SGI  workstation.  The  control  file,  which  includes  such  parameters  as  the  point  spread  function, 
bounding  box,  and  data  filename,  was  modified  to  test  the  algorithm  over  the  ultrasound 
images. 

4.5.2  Results  of  Genetic  Algorithm 

The  genetic  algorithm  technique  presented  by  Heckman  and  Krumm  was  tested  on 
several  ultrasound  images.  The  control  file  for  each  image  was  modified  to  use  the  correct 
bounding  boxes  for  the  location  of  the  bone  structures.  It  took  approximately  20  minutes  on  a 
SGI  Indigo  2  to  find  the  three  best-fitting  contours  for  the  outer  skin,  tibia,  and  fibula. 
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Figure  15  Contours  from  GA 


The  composite  image  of  the  resulting  contours  and  original  image  shown  in  Figure  15 
does  not  exhibit  smooth  boundary  contours.  The  presence  of  noise  within  the  image  and  the 
wide  intensity  variations  were  the  factors  that  caused  the  irregular  contours.  To  compensate,  the 
search  lines  used  to  align  the  contours  of  the  closed  cubic  spline  to  an  image's  intensity  peaks 
were  shortened.  With  this  modification,  the  GA  generated  smoother  contours  that  were  less 
sensitive  to  noise  and  clutter  in  the  images.  Figure  16  shows  a  new  GA  composite  image  with 
smoother  contours. 
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Figure  16  Smoother  Contours  with  Modified  GA  Parameter 


The  bounding  box  parameters  for  each  slice  were  manually  located  and  placed  in  a 
control  file  since  each  image  slice  varies  slightly.  This  step  was  later  automated  by  executing 
the  enhanced  2D  Bayesian  technique  on  each  image  slice  to  extract  bone  regions.  Then  the 
bounding  boxes  associated  with  the  outer  skin  and  each  bone  region  were  written  to  a  control 
file  for  use  by  the  GA  to  find  the  best-fitting  contours. 

4.6  Summary 

Unless  otherwise  noted,  all  source  code  was  implemented  with  the  C++  object-oriented 
programming  language  and  compiled  using  the  CC  Silicon  Graphics  (SGI)  compiler.  All  code 
was  compiled  and  tested  on  a  SGI  Indigo  2  workstation  equipped  with  a  250  MHz 
microprocessor  and  192  MB  of  RAM  running  the  IRIS  6.2  operating  system. 

All  segmentation  techniques  were  evaluated  using  8  to  10  ultrasound  images  from  the 
dataset  of  20  images.  This  subset  of  images  was  adequate  due  to  lack  of  complex  bone 


structures  and  similar  anatomical  characteristics  between  ultrasound  slices  of  a  BK  amputee's 
residual  limb.  Figure  17  gives  timing  results  for  the  segmentation  algorithms  presented  in  this 
chapter. 


Segmentation  Technique 

Average  Time  to  Process  an  Ultrasound  Slice 

Original  2D  Bayesian 

2.9  seconds 

Enhanced  2D  Bayesian 

15.6  seconds 

3D  Bayesian 

28.9  seconds 

Marr-Hildreth  and  Morphological  Operators 

8.6  seconds 

Genetic  Algorithm 

18.2  minutes 

Figure  17  Timing  Results  for  Segmentation 


The  four  major  segmentation  techniques  implemented  in  this  research  were  chosen 
based  on  their  appropriateness  and  potential  for  accurately  segmenting  ultrasound  images  of 
the  lower  limb.  Both  the  2D  multiresolution  and  3D  Bayesian  techniques  rely  on  a  region- 
based  segmentation  approach.  Both  techniques  use  an  adaptive  clustering  algorithm  that  takes 
into  account  spatial  constraints  and  local  intensity  changes.  The  2D  Bayesian  approach  required 
further  enhancement  to  accurately  distinguish  bone  region  in  the  ultrasound  images.  Like  the 
2D  Bayesian,  the  filtered  3D  approach  produced  similar  segmented  images  that  required  more 
postprocessing.  The  last  two  techniques  are  boundary-based  segmentation  approaches.  One 
combined  an  edge  detection  operator  with  morphological  operators.  The  other  used  a  model- 
driven  genetic  algorithm  that  searches  for  the  best-fitting  contours  using  a  fitness  function.  Out 
of  the  four  techniques,  the  enhanced  2D  Bayesian  produced  the  best  overall  segmentation  in 
terms  of  accurately  extracting  bone  structures  and  processing  speed.  The  two  techniques  with 
faster  processing  times  did  not  accurately  segment  the  bone  regions  in  the  ultrasound  images. 
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5  Generation  of  Surface  Polygonal  Database 


5.1  Introduction 

This  chapter  presents  a  method  using  public-domain  software  to  create  a  3D  lower  limb 
model  that  can  be  used  by  a  CASD  system  to  obtain  underlying  bone  structures.  Since  the 
enhanced  2D  Bayesian  provides  the  best  overall  segmentation  in  terms  of  accuracy  and 
processing  time,  this  was  the  only  technique  used  to  segment  and  create  a  realistic  lower  limb 
model.  After  segmentation  to  identify  boundaries  of  different  tissues  (the  skin  and  two  bones), 
the  Visualization  Toolkit  version  1.1  was  used  to  produce  the  surface  database  of  polygons 
(triangles)  and  render  the  lower  limb  model. 

The  Visualization  Toolkit  (VTK)  is  an  object-oriented  C++  library  of  data  objects  and 
algorithms  for  3D  rendering  written  by  William  Schroeder,  Bill  Lorenson  and  Ken  Martin 
[SCHR96].  The  authors  designed  the  library  specifically  for  scientific  and  medical  visualization. 
Their  goal  was  to  develop  the  library  for  easy  use  by  novices  and  experts  alike.  VTK  is  public 
domain  software  available  on  the  World  Wide  Web  or  a  CD  that  comes  with  the  book  written 
by  the  authors. 

VTK  uses  a  pipeline  architecture  to  process  and  manipulate  the  data  for  rendering 
[SCHR96].  Figure  18  illustrates  the  data  flow  of  the  architecture.  A  data  source  provides  the 
input  for  visualization  that  is  stored  as  a  VTK  data  object.  The  VTK  object  can  then  be  processed 
with  one  or  more  filters  such  as  "marching  cubes",  decimation,  or  contouring  filters.  The  data 
output  from  the  filter  is  also  stored  as  a  VTK  data  object.  The  data  can  be  mapped  to  geometric 
primitives  for  graphical  display  using  OpenGL  or  GL  on  the  SGI  workstations  or  written  out  to 
a  file. 
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Figure  18  Data  Flow  Through  VTK  Pipeline  Architecture 


5.2  3D  Surface  Rendering  Using  VTK 

For  this  research,  VTK's  "marching  cubes"  filter  was  chosen  to  create  a  lower  limb 
model  from  the  ultrasound  images.  The  "marching  cubes"  filter  identifies  the  polygons  cutting 
through  each  voxel  for  given  isovalues  of  each  tissue  class  and  creates  a  skin  contour  of  the 
original  image  dataset.  However,  accurate  extraction  of  bone  structures  requires  use  of 
segmented  image  data.  The  segmented  image  slices  are  better  for  creating  the  surface  contours 
since  there  are  continuous  boundary  definitions  and  constant  intensity  values  associated  with 
the  outer  skin  and  bone  structures. 

Before  VTK  could  be  used  to  decimate,  or  reduce  the  number  of  polygons  defining 
contours,  the  data  for  the  skin  and  bone  contours  has  to  be  processed  using  a  cleaning  filter  to 
eliminate  redundant  polygons  generated  by  the  "marching  cubes"  filter.  Then  to  eliminate 
noise  and  artifacts,  a  connectivity  filter  extracts  the  largest  connected  contour  belonging  to  the 
outer  skin  of  the  limb.  Afterwards,  the  data  is  piped  to  a  VTK  filter  for  checking  polygon 
normals.  This  process  ensures  that  each  polygon's  normal  correctly  identifies  the  front  and  back 
faces  of  the  skin  and  bone  surfaces.  Finally,  the  surface  model  of  the  limb  is  colored  and 
illuminated  with  VTK's  built-in  OpenGL  methods  to  create  a  realistic  3D  rendering  that  can  be 
rotated  and  zoomed  in  and  out  for  better  viewing.  Each  set  of  polygons  for  the  outer  skin  and 
bones  is  also  written  to  a  file  in  a  MOVIE.BYU  formate  to  be  imported  by  CARD  Laboratory's 
system.  Figure  19  diagrams  each  of  the  VTK  classes  created  to  generate  graphics  and  output 
data  to  MOVIE.BYU  file  format  of  the  lower  limb  model. 
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Figure  19  VTK  Classes  Called  to  Model  Lower  Limb 


An  additional  VTK  class  was  written  to  read  in  the  raw  binary  image  files  and  store  the 
data  as  a  VTK  data  structure.  A  special  VTK  class  was  also  written  to  output  the  results  of  the 
connectivity  filter  to  the  MOVIE.BYU  format.  The  executable  relies  on  the  use  of  a  configuration 
file  that  has  the  directory  path  of  the  data  files,  image  slice  range,  isovalue,  and  aspect  ratios  for 
the  X,  y,  and  z  dimensions. 
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53  Results  of  Using  VTK 


Usually,  a  3D  model  can  be  generated  using  the  "marching  cubes"  algorithm  on  the 
original  ultrasound  dataset.  However,  due  to  noise  and  widely  varying  intensities  throughout 
the  image  (even  within  similar  tissue  classes),  the  "marching  cubes"  technique  without  any 
preprocessing  did  not  result  in  an  accurate  3D  model.  Figure  20  shows  the  3D  representation 
obscured  by  noise  and  artifacts. 

The  model  produced  by  merely  using  the  segmented  image  slices  for  input  into  the 
VTK  procedure  did  not  produce  smooth  contours  for  the  outer  skin.  Due  to  step-like  intensity 
changes  at  the  boundaries,  the  skin  surface  contour  had  a  very  jagged  appearance.  Figure  21 
displays  the  lower  limb  model  from  segmented  images. 


Figure  20  3D  Model  Created  Using  Only  Marching  Cubes 
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Figure  21  3D  Model  of  Lower  Limb  from  Segmented  Images 


To  address  these  problems,  a  subsampled  set  of  images  was  created  using  MATLAB's 
imresize,  a  function  that  reduces  and  filters  an  image,  on  the  original  ultrasound  dataset.  This 
new  set  of  images  was  substituted  for  the  segmented  images.  Applying  all  the  VTK  classes 
shown  in  Figure  19  on  the  subsampled  images  created  polygons  defining  the  outer  skin 
contour.  The  subsampled  64  x  64  pixel  images  with  aspect  ratio  of  2.8  x  2.8  x  6  produced  more 
realistic,  smoother  skin  contours. 

Since  the  bones  structures  could  not  easily  be  detected  using  the  original  dataset,  the 
segmented  256  x  256  pixel  images  were  used  to  more  accurately  reconstruct  the  bone  structures 
with  aspect  ratio  of  0.7  x  0.7  x  6.  Although  this  produced  a  smooth  outer  skin  contour  and 
accurately  portrayed  underlying  bone  structures  (Figure  22  shows  a  new  lower  limb  model), 
there  were  several  tiny  holes  present  in  the  model.  These  holes  resulted  from  the  lack  cf 
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constant  intensity  values  for  the  boundary  definitions  in  the  images.  Since  the  ultrasound 
images  originally  had  wide  intensity  variations  with  bright  streaks  marking  the  boundaries, 
the  gaps  in  the  boundaries  produced  the  holes  present  in  the  model.  These  holes  did  not  result 
with  use  of  segmented  images  since  they  have  specific  intensity  values  associated  with  different 
tissue  regions.  Therefore,  VTK's  "marching  cubes"  filter  a  connected  set  of  polygons  that 
completely  define  the  skin  contour. 

5.4  Summary 

The  VTK  library  created  by  Schroeder,  Lorensen  and  Martin  provides  a  fast  object- 
oriented  architecture  for  visualizing  3D  medical  images.  Various  processing  filters  from  the 
VTK  library  created  a  3D  lower  limb  model  from  the  ultrasound  image  dataset.  However,  the 
limitations  presented  by  the  quality  of  the  ultrasound  imaging  technique  contributed  to  initial 
problems  in  creating  a  3D  representation  of  a  amputee's  residual  limb.  To  create  a  lower  limb 
model  that  clearly  showed  the  skin  and  bone  contours,  VTK  used  the  segmented  dataset  to 
generate  the  surface  polygons.  The  enhanced  2D  Bayesian  algorithm  was  selected  to  segment 
the  entire  ultrasound  dataset  based  on  its  performance  as  discussed  in  the  previous  chapter. 
Then  the  filter  classes  diagramed  in  Figure  19  were  used  to  generate  a  realistic  3D  lower  limb 
model  of  the  outer  skin  and  bone  contours.  After  creating  a  realistic  3D  lower  limb  model,  the 
surface  polygons  were  written  to  a  file  in  a  MOVIE.BYU  format.  Any  application  supporting 
MOVIE.BYU  files  can  import  the  3D  lower  limb  model. 
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Figure  22  3D  Model  of  Lower  Limb  with  Smooth  Outer  Contour 
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6  Conclusion  and  Recommendations 


Ultrasound  imaging  is  quickly  gaining  popularity  because  it  is  efficient,  inexpensive, 
and  non-invasive  [NELS93].  However,  to  produce  rendered  images  that  have  the  high  quality 
of  other  modalities,  more  research  is  needed  to  develop  faster,  more  accurate  techniques  that 
meet  the  special  challenges  poised  by  the  unique  low-noise  ultrasound  images  described  in 
Chapter  2.  Even  with  image  processing  techniques,  problems  associated  with  noise,  incomplete 
boundaries,  and  wide  intensity  variations  within  images  still  present  obstacles  for  segmentation 
and  creation  of  an  accurate  3D  lower  limb  model.  These  limitations  are  never  completely 
eliminated  during  image  processing  since  the  goal  is  to  optimize  image  quality  while 
preserving  important  details. 

Segmentation  is  very  important  for  creating  accurate  models  for  use  with  a  CASD 
system.  My  research  identified  and  evaluated  four  potential  segmentation  techniques  and 
enhancements  that  can  be  used  to  segment  low-noise  ultrasound  residual  limb  images.  I 
developed  object-oriented  C++  code  or  MATLAB  code  to  test  each  segmentation  algorithm.  I 
also  used  the  Visualization  Toolkit  to  generate  a  realistic  3D  model  based  on  the  segmented 
image  dataset  for  use  in  a  computer-aided  socket  design  system. 

Although  this  research  was  partially  focused  on  segmenting  the  outer  skin  contour, 
which  was  done  without  difficulty,  the  important  research  dealt  with  the  extraction  of  the  bone 
structures.  The  enhanced  2D  Bayesian  algorithm  produced  fairly  accurate  segmentation  of  the 
lower  limb.  Dr.  Ping  He  and  Dr.  Kefu  Xue,  two  biomedical  experts  from  Wright  State 
University,  verified  the  accuracy  of  segmentation  of  the  tibia,  the  larger  bone.  Dr.  He  stated 
that  for  prosthetic  socket  design,  the  location  of  the  tibia  is  the  most  important  since  it  tends  to 
bear  most  of  the  weight  and  requires  relief  areas.  The  fibula,  the  smaller  bone,  was  much 
harder  to  distinguish  due  to  the  lack  of  boundary  definitions  from  the  tibia's  backshadow.  The 
fibula  also  needs  accurate  segmentation  because  prosthetists  may  need  to  know  of  its  interaction 
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with  the  skin  and  tibia.  Therefore,  more  research  is  need  to  overcome  the  difficulties  inherent 
in  ultrasound  imaging  to  accurately  extract  the  smaller  bone  region.  The  3D  Bayesian 
algorithm  did  not  produce  more  accurate  results  than  the  enhanced  2D  Bayesian  and  processed 
each  slice  slower  than  the  2D  approach.  Although  the  Marr-Hildreth  algorithm  was  slightly 
faster,  the  its  accuracy  was  also  less  than  that  of  the  2D  Bayesian. 

Segmentation  may  be  improved  by  creating  a  model-driven  approach,  such  as  the 
genetic  algorithm  technique  presented  by  Heckman  and  Krumm.  This  type  of  approach  could 
help  offset  the  distractions  caused  by  noise  since  their  GA  was  able  to  produce  accurate 
segmentation  of  both  regions.  However,  Heckman  and  Krumm's  technique  is  not  automated 
and  relies  on  some  human  interaction  for  locating  bounding  boxes  of  the  bone  structures. 
Although  this  step  could  be  automated  using  the  2D  Bayesian  algorithm  to  extract  the  bounding 
boxes,  this  additional  step  would  increase  the  overhead.  The  genetic  algorithm  is  quite  slow 
compared  to  the  other  algorithms  since  it  covers  a  huge  search  space  and  adding  the  extra  step 
to  automate  segmentation  would  further  decrease  speed.  Furthermore,  each  execution  of  the 
algorithm  produces  different  results  since  the  genetic  algorithm  is  based  on  finding  an  optimal 
solution  not  the  best  one. 

Overall,  for  quick  speed  and  relative  accuracy  the  enhanced  2D  Bayesian  produced  the 
best  segmentation  of  ultrasound  images.  However,  further  work  is  required  to  increase 
accuracy  in  the  segmentation  of  the  two  bones  especially  the  fibula.  Future  success  may  result 
from  trying  to  incorporate  and  mimic  the  human  perceptucJ  system  and  expert  knowledge  of 
ultrasound  specialist  and  prosthetists. 

Since  the  scope  of  this  research  focused  primarily  on  segmentation,  the  improvement  of 
surface  tracking  techniques  was  not  researched  in  great  detail.  However,  if  increased  speed  is 
important,  faster  surface  tracking  could  result  by  modifying  the  "marching  cubes"  algorithm  in 
the  VTK  library.  The  algorithm  could  be  enhanced  by  using  the  approach  proposed  by  Livnat, 
et  al  Unlike  the  "marching  cubes"  algorithm  which  inspects  each  and  every  voxel  to  see  if  a 
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surface  intersects  it,  their  near-optimal  isosurface  extraction  (NOISE)  algorithm,  a  faster 
technique  released  in  March  1996,  uses  a  multidimensional  binary  tree  structure  to  quickly 
eliminate  voxels  that  do  not  meet  the  criteria  for  being  intersected.  As  a  result,  the  algorithm  is 
faster  and  more  efficient  in  identifying  the  voxels  that  have  surfaces  cutting  through  them 
[LIVN96]. 

In  addition,  the  primary  concern  of  my  research  was  on  the  evaluation  of  various 
segmentation  techniques  on  specialized  ultrasound  images.  Therefore,  I  mainly  concentrated 
on  fast  implementation  of  the  different  algorithms  and  did  not  emphasize  speed  and  efficiency 
of  processing.  Therefore,  my  implementation  of  the  different  segmentation  techniques  may 
need  further  analysis  for  modifications  leading  to  faster  processing  and  lower  memory 
requirements. 

With  further  research  and  improvements  on  the  work  presented  here,  a  more  realistic 
and  accurate  model  of  the  residual  limb  of  a  BK  amputee  could  be  created  to  help  solve  the 
problems  that  exist  with  conventional  prosthetic  socket  design  which  would  save  money, 
provide  greater  comfort,  and  more  mobility  to  amputees. 
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Appendix  A:  Original  Ultrasound  Images  Slices 
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Ultrasound  Image  Slices  1  -  4  of  3D  Dataset 


Ultrasound  Image  Slices  5  -  8  of  3D  Dataset 
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Ultrasound  Image  Slices  9  - 12  of  3D  Dataset 


57 


Ultrasound  Image  Slices  13  - 16  of  3D  Dataset 
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Ultrasound  Image  Slices  17  -  20  of  3D  Dataset 
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Appendix  B:  2D  Bayesian  Segmented  Images 
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Ultrasound  Segmented  Slices  5  -  8  of  3D  Dataset 
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Ultrasound  Segmented  Slices  9  - 12  of  3D  Dataset 
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Ultrasound  Segmented  Slices  13  - 16  of  3D  Dataset 
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Appendix  C:  Enhanced  2D  Bayesian  Segmented 

Images 
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Ultrasound  Segmented  Slices  17  -  20  of  3D  Dataset 
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